An Ionosphere Estimation Algorithm for WAAS Based on Kriging

GPS alone cannot provide the integrity needed for air navigation. Several error sources deteriorate the precision of the position estimate. One of the largest and more unpredictable sources of error for single frequency users is the ionosphere. For this reason, ionospheric behavior drives the performance of the Wide Area Augmentation System (WAAS). At any given time, the only information we have of the ionosphere is a limited amount of Total Electron Content (TEC) measurements. As a consequence, in order to estimate the ionospheric delay and get a confidence bound on such an estimate, we need to understand the spatial structure of the ionosphere over the region of interest. Using the thin shell model framework, where each TEC measurement is identified as a location on the thin shell, labeled the Ionospheric Pierce Point (IPP), the problem is reduced to a 2-dimensional problem. Once we have a good description of a nominal ionosphere, there are two questions that need to be answered before estimating the delay at a given IPP: Are the IPP measurements compatible with the assumed nominal model of the ionosphere? How relevant are the IPP measurements to the location we need to estimate? To answer the first one, an accurate characterization of the ionosphere in nominal conditions is needed. The large observed stationarity violations make this latter question very difficult. A worst case based method to determine the spatial structure of the nominal ionosphere in terms of the variogram, or, equivalently, the covariance is presented. The technique called ‘kriging’ produces at each location an estimate and a confidence bound on the estimate, the kriging variance. The particular behavior of the kriging variance at the edge of coverage allows us to intuitively define the ‘well sampled’ region. We show that a carefully designed estimation algorithm based on kriging could provide confidence bounds on the ionospheric delay corrections allowing WAAS to meet the GNSS Landing System requirements. INTRODUCTION WAAS corrects, among other sources of error, for the user’s ionospheric delay errors and places strict confidences on those corrections under all conditions [1]. The estimation process is simplified by making the thin shell approximation [2]. Each ray path has a corresponding Ionospheric Pierce Point (IPP), where that path intersects the thin shell. Each ionospheric TEC measurement is transformed in a Vertical Ionospheric Delay. The WAAS correction message specifies the vertical delay values as well as the confidence at each node, labeled σGIVE [3]. The user interpolates these values to find the ionospheric delays corresponding to its satellites in view. The current algorithm is described in [4]. In order to increase the performance of the system, in particular to be able to meet GNSS Landing System requirements, the confidence bound on the correction would need to be reduced by more than a factor of two: from a σGIVE of 1.5 m to .6 m. In fact, under nominal conditions, any reasonable estimation algorithm yields excellent estimates. Unfortunately, the existence of sudden changes in the ionosphere, from quiet to stormy conditions [5], [6], [7], together with the fact that the ionosphere is irregularly sampled, has forced the system to be overly conservative most of the time. A storm detector based on the chisquare test, and a metric measuring the degree of coverage take care of the irregularity threats [4], [8]. It is mainly in the characterization of the undersampled regions that there is room for improvement. It was shown in [9] that the technique called kriging –a particular minimum mean-square estimator – could mitigate the undersampled problem. Kriging has long been used for mapping purposes, because it provides a natural interpolation. As such it has been applied to map TEC measurements as well as other ionospheric parameters [10], [11]. Kriging has also many attractive features for WAAS [9], for both the definition of the ‘well sampled’ region and the confidence bound generation. However, a critical step in the estimation when using kriging is the knowledge of the underlying deterministic and random spatial structure [12]. In this work, we are interested in the limitations of ionospheric delay estimation due to ionospheric behavior and due to the location of the IPP measurements. We will therefore omit the difficulties introduced by the bandwidth limitation [4] by assuming that the user has the same information as the master station, that is, all the IPP locations and measurements. For the validation of the proposed solutions, we will also consider that the measurements have very low noise, as it is the case in WAAS ‘supertruth’ data. In this set of data, the redundancy of receivers and the post-processing allows the isolation and removal of faults [5]. However, it will be explained how to include the higher level of noise that real time measurements have. In the first part, we will show a method to refine the spatial characterization of the ionosphere in nominal conditions, and come up with a conservative model. In the second part we will describe the algorithm that is optimal under the characterization given in the first part; we will also show how to heuristically generate the ‘well sampled’ region from the kriging variance. We finally will present the cross validation results over several days worth of supertruth data, which includes the most extreme ionospheric conditions observed during the current solar cycle [7]. SPATIAL STRUCTURE OF THE VERTICAL IONOSPHERIC DELAY A useful and intuitive way of characterizing the spatial structure of a random field is the variogram [13]. The variogram measures the loss of confidence of a random function as we depart from the measurements, assuming that there is no deterministic underlying trend, and that the random function is stationary. It is defined as: ( ) ( ) ( ) ( ) { } 2 12 h E I u h I u γ = + − (1) where u is a location, h is a vector and I is the random function that we want to estimate; in our case, it is the vertical ionospheric delay. At each epoch the variogram was estimated by computing all the difference between pairs of measurements Iv(ui)-Iv(uj) and the corresponding distance |ui-uj| [5], [9]. The pairs were classified according to their distance in different categories. In this study we used the bins [0 300 km], [300 km 600 km] etc. For each bin we estimated the variance of the resulting empirical distribution. Since we wanted to get a conservative estimate of the variance, we used the gaussian overbound [14] of the empirical distribution. Figure 1. Variograms for a quiet day Figure 1. shows several experimental variograms for a quiet day. The flat variograms correspond to the ionosphere during the night, and the more variable correspond to daytime. In fact, the ionospheric behavior is well approximated by a planar trend [5]. It is the planar trend that is responsible for most of the variability, so it does not correspond to random variability. We can have a better idea of the range of the trend by examining a variogram obtained through the classical formula –instead of computing the gaussian overbound in each bin, we compute the mean of the squares[9], [12] which is smoother than the variogram obtained through the gaussian overbound. Figure 2 shows such a variogram. Figure 2. Variogram for a quiet day obtained with the classical formula. Although we cannot use this variogram to bound the ionospheric behavior, we can extract some qualitative information. The overall shape is determined by the presence of a trend, which produces the parabolic behavior from 0 to 2000km. However, we see that for lags larger than 2500 km, the planar trend is not as clear. The other very important characteristic of the variogram is the behavior at the origin [12]. In this variogram, at small lags we have a non zero derivative. This means that even if we remove the planar trend there is random spatial structure left [12], [13]. We can express this idea by writing that, locally, the vertical ionosphere delay is given by: ( ) 0 1 2 , ( , ) v I x y a a x a y R x y = + + + (2) where the first three terms define the planar trend and R(x,y) is a random function with zero mean and a given variogram. Now, we would like to characterize the random behavior of R, that is, to find a variogram that describes it conservatively. De-trending the data is a classical problem in spatial statistics [15]. It is problematic to just fit a plane to the data and compute the variogram of the residuals. Some measurements might have high leverage on the fit, either due to their large value or their geometry [16], contaminating the original data. For example a measurement that represents an outlier to the plane will be smoothed if we include it in the planar fit; in our description of the ionosphere, we want to be aware if such outliers do exist, and we want to understand their spatial dependency. We therefore need a method that does not contaminate the original data, or that contaminates it in a conservative way. The common practice in spatial statistics is to compare only the pairs of measurements that lie in a direction unaffected by the trend [12], [15]. This way we ensure that the original data is not affected by the fit. For our purposes, there are two problems in this approach. First, we end up with a very small amount of measurements, thus getting less statistical significance. Second, we assume that the variability in the direction of the trend is the same as in the direction orthogonal to the trend. These two problems go against our worst case approach. Instead, we decided to proceed in the following way: for each pair of measurements, we did a planar fit on the remaining measurements (up to a certain radius) and assumed this planar fit would describe the trend. The advantage of this approach is that the data used to generate the pairs will not influence the trend. Also, we do not limit the number of pairs. The inconvenience is that we do not take into account the uncertainty that the plane might have, due to the geometry of the measured locations, thus not defining correctly the trend. However, this is likely to bias the variogram in a conservative way. We mitigated this effect by considering only the densest regions in terms of IPP measurements. Figure 3. Experimental residual variograms for a quiet day. Figure 3 shows a series of experimental residual variograms generated using this method. We used only pairs of measurements with latitudes 30 deg and 45 deg, and longitudes between -105 deg and -85 deg. The radius of search used was 1500 km. The nominal day used was July 2, 2000. If the random function R had no spatial dependency, we would expect those variograms to be flat. We see instead that there is a clear spatial structure once the trend is removed. Since we are going to be using the variogram for estimation, we need to find an admissible variogram [11], [13] that is close enough to the experimental variogram. The simplest we can find is a linear model with a non zero intercept at the origin. However, such a model cannot be described by a covariance, because the variogram has no finite limit as the separation becomes infinite [12] –we will see later that we need the covariance. We chose instead a model variogram with a linear behavior at the origin and a finite limit at infinity –a sill. An admissible variogram having these properties is given by: